#include<iostream>
#include<vector>
#include<queue>
#include<malloc.h>
#include <math.h>
#include <lapacke.h>
#include <algorithm>
#include <iterator>



void BDFs(int order)

{ 
  int i,j,n,q;
  double u,k1,k2,T1,T2,a4,a3,a2,a1,a0,b4,y1,y2,tol;
  double c11,c12,c13,c14,c21,c22,c23,c24,c31,c32,c33,c34,c41,c42,c43,c44,d1,d2,d3,d4;
  double *u1,*u2,*u4,*u5,v1,v2,v4,v5;
  double e1,e2,emaxnorm,premaxnorm,Crate;
  double start,finish,duration;

  u1=(double *)malloc(sizeof(double)*10000000);
  u2=(double *)malloc(sizeof(double)*10000000);
  u4=(double *)malloc(sizeof(double)*10000000);
  u5=(double *)malloc(sizeof(double)*10000000);

  u1[0]=0.994;u2[0]=0;u4[0]=0;u5[0]=-2.0015851063790825224;u=1/81.45;
  T1=17.06521656015796;T2=19.14045706162071;
  k1=T1/n;

  if (order==1){
    a4=1.0;a3=-1.0;a2=0;a1=0;a0=0;b4=1.0;}
  else if (order==2){
    a4=1.0;a3=-4.0/3.0;a2=1.0/3.0;a1=0;a0=0;b4=2.0/3.0;}
  else if (order==3){
    a4=1.0;a3=-18.0/11.0;a2=9.0/11.0;a1=-2.0/11.0;a0=0;b4=6.0/11.0;}
  else{
    a4=1.0;a3=-48.0/25.0;a2=36.0/25.0;a1=-16.0/25.0;a0=3.0/25.0;b4=12.0/25.0;}

  std::cout<<"p="<<order<<std::endl;
  std::cout<<"the BDFs method(situation1):"<<std::endl;
      std::cout<<"  节点数  "<<"  运行时间  "<<" 误差的最大范数  "<<"  收敛阶  "<<std::endl;
    
  //
 
  
  n=5e3;k1=T1/n;
  for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }

  for (i=4;i<=n;i++){
    v1=u1[i-1];
    v2=u2[i-1];
    v4=u4[i-1];
    v5=u5[i-1];
    
    //求右侧向量;
    d1=a0*u1[i-4]+a1*u1[i-3]+a2*u1[i-2]+a3*u1[i-1]+a4*u1[i-1]-k1*b4*u4[i-1];
    d2=a0*u2[i-4]+a1*u2[i-3]+a2*u2[i-2]+a3*u2[i-1]+a4*u2[i-1]-k1*b4*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    d3=a0*u4[i-4]+a1*u4[i-3]+a2*u4[i-2]+a3*u4[i-1]+a4*u4[i-1]-k1*b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    d4=a0*u5[i-4]+a1*u5[i-3]+a2*u5[i-2]+a3*u5[i-1]+a4*u5[i-1]-k1*b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
    c11=a4;c12=0;c13=-k1*b4;c14=0;
    c21=0;c22=a4;c23=0;c24=-k1*b4;
    y1=pow(y1,5.0/3.0);y2=pow(y2,5.0/3.0);
    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);c33=a4;c34=-2*k1*b4;
    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);
    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c43=2*k1*b4;c44=a4;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
    for (j=1;j<=100;j++){
      if (tol>1e-6){
	u1[i-1]=u1[i];
	u2[i-1]=u2[i];
	u4[i-1]=u4[i];
	u5[i-1]=u5[i];
        //newton迭代
	//求右侧向量;
    d1=a0*u1[i-4]+a1*u1[i-3]+a2*u1[i-2]+a3*v1+a4*u1[i-1]-k1*b4*u4[i-1];
    d2=a0*u2[i-4]+a1*u2[i-3]+a2*u2[i-2]+a3*v2+a4*u2[i-1]-k1*b4*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    d3=a0*u4[i-4]+a1*u4[i-3]+a2*u4[i-2]+a3*v4+a4*u4[i-1]-k1*b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    d4=a0*u5[i-4]+a1*u5[i-3]+a2*u5[i-2]+a3*v5+a4*u5[i-1]-k1*b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
    c11=a4;c12=0;c13=-k1*b4;c14=0;
    c21=0;c22=a4;c23=0;c24=-k1*b4;
    y1=pow(y1,5.0/3.0);y2=pow(y2,5.0/3.0);
    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);c33=a4;c34=-2*k1*b4;
    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);
    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c43=2*k1*b4;c44=a4;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
      }
      else{
	break;}
    }
    u1[i-1]=v1;
    u2[i-1]=v2;
    u4[i-1]=v4;
    u5[i-1]=v5;
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);

  //
  for (q=1;q<=10;q++){
  premaxnorm=emaxnorm;
  start = clock();    
  n=n*2;
  k1=T1/n;
  
  for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }

  for (i=4;i<=n;i++){
    v1=u1[i-1];
    v2=u2[i-1];
    v4=u4[i-1];
    v5=u5[i-1];
    
    //求右侧向量;
    d1=a0*u1[i-4]+a1*u1[i-3]+a2*u1[i-2]+a3*u1[i-1]+a4*u1[i-1]-k1*b4*u4[i-1];
    d2=a0*u2[i-4]+a1*u2[i-3]+a2*u2[i-2]+a3*u2[i-1]+a4*u2[i-1]-k1*b4*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    d3=a0*u4[i-4]+a1*u4[i-3]+a2*u4[i-2]+a3*u4[i-1]+a4*u4[i-1]-k1*b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    d4=a0*u5[i-4]+a1*u5[i-3]+a2*u5[i-2]+a3*u5[i-1]+a4*u5[i-1]-k1*b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
    c11=a4;c12=0;c13=-k1*b4;c14=0;
    c21=0;c22=a4;c23=0;c24=-k1*b4;
    y1=pow(y1,5.0/3.0);y2=pow(y2,5.0/3.0);
    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);c33=a4;c34=-2*k1*b4;
    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);
    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c43=2*k1*b4;c44=a4;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
    for (j=1;j<=100;j++){
      if (tol>1e-6){
	u1[i-1]=u1[i];
	u2[i-1]=u2[i];
	u4[i-1]=u4[i];
	u5[i-1]=u5[i];
        //newton迭代
	//求右侧向量;
    d1=a0*u1[i-4]+a1*u1[i-3]+a2*u1[i-2]+a3*v1+a4*u1[i-1]-k1*b4*u4[i-1];
    d2=a0*u2[i-4]+a1*u2[i-3]+a2*u2[i-2]+a3*v2+a4*u2[i-1]-k1*b4*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    d3=a0*u4[i-4]+a1*u4[i-3]+a2*u4[i-2]+a3*v4+a4*u4[i-1]-k1*b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    d4=a0*u5[i-4]+a1*u5[i-3]+a2*u5[i-2]+a3*v5+a4*u5[i-1]-k1*b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
    c11=a4;c12=0;c13=-k1*b4;c14=0;
    c21=0;c22=a4;c23=0;c24=-k1*b4;
    y1=pow(y1,5.0/3.0);y2=pow(y2,5.0/3.0);
    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);c33=a4;c34=-2*k1*b4;
    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);
    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c43=2*k1*b4;c44=a4;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
      }
      else{
	break;}
    }
    u1[i-1]=v1;
    u2[i-1]=v2;
    u4[i-1]=v4;
    u5[i-1]=v5;
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);
  Crate=premaxnorm/emaxnorm;Crate=log2(Crate);
  finish = clock();
  duration = (double)(finish - start) / CLOCKS_PER_SEC;
  std::cout<<"  "<<n<<"     " <<duration<<"    "  <<emaxnorm<<"      "<<Crate  <<std::endl;
 
  }

 T1=T2,u1[0]=0.87978;u2[0]=0;u4[0]=0;u5[0]=-0.3797;
      std::cout<<"the BDFs method(situation2):"<<std::endl;
      std::cout<<"  节点数  "<<"  运行时间  "<<" 误差的最大范数  "<<"  收敛阶  "<<std::endl;

n=1250;k1=T1/n;
  for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }

  for (i=4;i<=n;i++){
    v1=u1[i-1];
    v2=u2[i-1];
    v4=u4[i-1];
    v5=u5[i-1];
    
    //求右侧向量;
    d1=a0*u1[i-4]+a1*u1[i-3]+a2*u1[i-2]+a3*u1[i-1]+a4*u1[i-1]-k1*b4*u4[i-1];
    d2=a0*u2[i-4]+a1*u2[i-3]+a2*u2[i-2]+a3*u2[i-1]+a4*u2[i-1]-k1*b4*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    d3=a0*u4[i-4]+a1*u4[i-3]+a2*u4[i-2]+a3*u4[i-1]+a4*u4[i-1]-k1*b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    d4=a0*u5[i-4]+a1*u5[i-3]+a2*u5[i-2]+a3*u5[i-1]+a4*u5[i-1]-k1*b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
    c11=a4;c12=0;c13=-k1*b4;c14=0;
    c21=0;c22=a4;c23=0;c24=-k1*b4;
    y1=pow(y1,5.0/3.0);y2=pow(y2,5.0/3.0);
    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);c33=a4;c34=-2*k1*b4;
    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);
    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c43=2*k1*b4;c44=a4;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
    for (j=1;j<=100;j++){
      if (tol>1e-6){
	u1[i-1]=u1[i];
	u2[i-1]=u2[i];
	u4[i-1]=u4[i];
	u5[i-1]=u5[i];
        //newton迭代
	//求右侧向量;
    d1=a0*u1[i-4]+a1*u1[i-3]+a2*u1[i-2]+a3*v1+a4*u1[i-1]-k1*b4*u4[i-1];
    d2=a0*u2[i-4]+a1*u2[i-3]+a2*u2[i-2]+a3*v2+a4*u2[i-1]-k1*b4*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    d3=a0*u4[i-4]+a1*u4[i-3]+a2*u4[i-2]+a3*v4+a4*u4[i-1]-k1*b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    d4=a0*u5[i-4]+a1*u5[i-3]+a2*u5[i-2]+a3*v5+a4*u5[i-1]-k1*b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
    c11=a4;c12=0;c13=-k1*b4;c14=0;
    c21=0;c22=a4;c23=0;c24=-k1*b4;
    y1=pow(y1,5.0/3.0);y2=pow(y2,5.0/3.0);
    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);c33=a4;c34=-2*k1*b4;
    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);
    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c43=2*k1*b4;c44=a4;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
      }
      else{
	break;}
    }
    u1[i-1]=v1;
    u2[i-1]=v2;
    u4[i-1]=v4;
    u5[i-1]=v5;
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);

  //
  for (q=1;q<=10;q++){
  premaxnorm=emaxnorm;
  start = clock();    
  n=n*2;
  k1=T1/n;
  
  for (i=1;i<=3;i++){
    u1[i]=u1[i-1]+k1*u4[i-1];
    u2[i]=u2[i-1]+k1*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    u4[i]=u4[i-1]+k1*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    u5[i]=u5[i-1]+k1*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
  }

  for (i=4;i<=n;i++){
    v1=u1[i-1];
    v2=u2[i-1];
    v4=u4[i-1];
    v5=u5[i-1];
    
    //求右侧向量;
    d1=a0*u1[i-4]+a1*u1[i-3]+a2*u1[i-2]+a3*u1[i-1]+a4*u1[i-1]-k1*b4*u4[i-1];
    d2=a0*u2[i-4]+a1*u2[i-3]+a2*u2[i-2]+a3*u2[i-1]+a4*u2[i-1]-k1*b4*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    d3=a0*u4[i-4]+a1*u4[i-3]+a2*u4[i-2]+a3*u4[i-1]+a4*u4[i-1]-k1*b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    d4=a0*u5[i-4]+a1*u5[i-3]+a2*u5[i-2]+a3*u5[i-1]+a4*u5[i-1]-k1*b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
    c11=a4;c12=0;c13=-k1*b4;c14=0;
    c21=0;c22=a4;c23=0;c24=-k1*b4;
    y1=pow(y1,5.0/3.0);y2=pow(y2,5.0/3.0);
    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);c33=a4;c34=-2*k1*b4;
    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);
    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c43=2*k1*b4;c44=a4;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
    for (j=1;j<=100;j++){
      if (tol>1e-6){
	u1[i-1]=u1[i];
	u2[i-1]=u2[i];
	u4[i-1]=u4[i];
	u5[i-1]=u5[i];
        //newton迭代
	//求右侧向量;
    d1=a0*u1[i-4]+a1*u1[i-3]+a2*u1[i-2]+a3*v1+a4*u1[i-1]-k1*b4*u4[i-1];
    d2=a0*u2[i-4]+a1*u2[i-3]+a2*u2[i-2]+a3*v2+a4*u2[i-1]-k1*b4*u5[i-1];
    y1=u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1);y1=pow(y1,1.5);
    y2=u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u);y2=pow(y2,1.5);
    d3=a0*u4[i-4]+a1*u4[i-3]+a2*u4[i-2]+a3*v4+a4*u4[i-1]-k1*b4*(2*u5[i-1]+u1[i-1]-(u*(u1[i-1]+u-1))/y1-((1-u)*(u1[i-1]+u))/y2);
    d4=a0*u5[i-4]+a1*u5[i-3]+a2*u5[i-2]+a3*v5+a4*u5[i-1]-k1*b4*(u2[i-1]-2*u4[i-1]-(u*u2[i-1])/y1-((1-u)*u2[i-1])/y2);
    d1=-d1;d2=-d2;d3=-d3;d4=-d4;

    //求雅克比矩阵;
    c11=a4;c12=0;c13=-k1*b4;c14=0;
    c21=0;c22=a4;c23=0;c24=-k1*b4;
    y1=pow(y1,5.0/3.0);y2=pow(y2,5.0/3.0);
    c31=-k1*b4*(1-u*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(u2[i-1]*u2[i-1]-2*(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c32=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);c33=a4;c34=-2*k1*b4;
    c41=-k1*b4*(3*u*u2[i-1]*(u1[i-1]+u-1)/y1+3*(1-u)*(u1[i-1]+u)*u2[i-1]/y2);
    c42=-k1*b4*(1-u*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u-1)*(u1[i-1]+u-1))/y1-(1-u)*(-2*u2[i-1]*u2[i-1]+(u1[i-1]+u)*(u1[i-1]+u))/y2);
    c43=2*k1*b4;c44=a4;

    //求解方程;
    double C[16]={c11,c21,c31,c41,
		  c12,c22,c32,c42,
		  c13,c23,c33,c43,
		  c14,c24,c34,c44};
    double D[4]={d1,d2,d3,d4};
    int ipiv[4];
    LAPACKE_dgesv(LAPACK_COL_MAJOR, 4, 1, C, 4, ipiv, D, 4);
    u1[i]=u1[i-1]+D[0];
    u2[i]=u2[i-1]+D[1];
    u4[i]=u4[i-1]+D[2];
    u5[i]=u5[i-1]+D[3];
    tol=fabs(D[0])+fabs(D[1])+fabs(D[2])+fabs(D[3]);
      }
      else{
	break;}
    }
    u1[i-1]=v1;
    u2[i-1]=v2;
    u4[i-1]=v4;
    u5[i-1]=v5;
  }
  e1=u1[n]-u1[0];
  e2=u2[n]-u2[0];
  e1=fabs(e1);e2=fabs(e2);
  emaxnorm=std::max(e1,e2);
  Crate=premaxnorm/emaxnorm;Crate=log2(Crate);
  finish = clock();
  duration = (double)(finish - start) / CLOCKS_PER_SEC;
  std::cout<<"  "<<n<<"     " <<duration<<"    "  <<emaxnorm<<"      "<<Crate  <<std::endl;
 
  }

  
}
	
	
      
    
    
    
  
